Introduction¶
Antimicrobial resistance is becoming a worldwide problem
Many pathogens are becoming resistant for many different types of antibiotics
The institute of advanced studies (IAS) is doing research to the spreading of gonorrhoea in Amsterdam
For gonorrhoeae there might be a new antibiotic available within years
It is thus important to find a suitable management plan to prevent new resistance
Mathematical model¶
- Aim of our research was to build a network model to model the disease spreading
- Our network models are based on an existing mathematical model by Fingerhurth et al. (2006) [1]
- Simple network model
- Stochastic network
Why use a network model?¶
- Sexual network portrays the sexual inter-relationships within a group of people and it is possible to study the spread of STIs in the context of social interactions [2]
- ODEs, in contrast, fail to capture the complexity of connections between cases
- By using a network a heterogeneous population can be introduced, which adds more detail in the properties of individuals. This can for example be useful in identifying high-risk groups.
- Adding this heterogeneity can have influence on the disease dynamics [3]
- Network properties can for example derived from data, in Amsterdam this is for example possible with GGD data
Mathematical model of antibiotic resistance N. gonorrhoeae spread¶
In [235]:
Image(filename= "math_model.png")
Out[235]:
In [242]:
plt.figure()
# plt.plot(t, sol_msm[1,:, 0], label='$S_{low}$')
# plt.plot(t, sol_msm[1,:, 1], label='$S_{high}$')
plt.plot(t, sol_hmw[1,:, 2], label='$I_{res, low}$')
plt.plot(t, sol_hmw[1,:, 3], label='$I_{res, high}$')
plt.plot(t, sol_hmw[1,:, 4], label='$I_{sen, low}$')
plt.plot(t, sol_hmw[1,:, 5], label='$I_{sen, high}$')
plt.legend(loc='best')
plt.xlabel('t (years)')
plt.title("Epidemic curve HMW group")
plt.ylabel("Fraction of population")
# plt.ylim(0, 1)
plt.xlim(0, 30)
plt.grid()
plt.show()
In [243]:
plt.figure()
ratio_sol_msm = (sol_msm[:,:, 2] + sol_msm[:,:, 3])/(sol_msm[:,:, 2] + sol_msm[:,:, 3] + sol_msm[:,:, 4] + sol_msm[:,:, 5])
plt.plot(t, ratio_sol_msm[1,:], label = "MSM", c = "g",linewidth=1)
plt.fill_between(t,ratio_sol_msm[0,:], ratio_sol_msm[2,:], alpha=0.5, facecolor = "g")
ratio_sol_hmw = (sol_hmw[:,:, 2] + sol_hmw[:,:, 3])/(sol_hmw[:,:, 2] + sol_hmw[:,:, 3] + sol_hmw[:,:, 4] + sol_hmw[:,:, 5])
plt.plot(t, ratio_sol_hmw[1,:], label = "HMW", c = "b",linewidth=1)
plt.fill_between(t,ratio_sol_hmw[0,:], ratio_sol_hmw[2,:], alpha=0.5, facecolor = "b")
plt.legend(loc='best')
plt.xlabel('t (years)')
plt.ylabel("Fraction of infected individuals")
plt.grid()
plt.ylim(0,1)
plt.show()
Simple Contact network¶
- We used a simplified version of the rules in the mathematical model
- The network is a random erdos renyi network with 500 nodes and an average degree of 10
The parameters
- Spontaneous healing rate $\nu$
- Treatment rate $\tau$
- Hosts develop resistance during treatment with probability $\mu$
- Transmision probability $\beta$
Build a network¶
In [264]:
# # create power law graph
# exp = 2.4
# seed = 123456789
# num = 1000
# graph = power_law_graph(exp, num, seed)
# or, create random network
degree = 10
num = 500
p = degree/num
seed = 1234
graph = random_network(p, num, seed)
# compute average degree
degrees = dict(graph.degree())
sum_of_edges = sum(degrees.values())
average_deg = sum_of_edges / len(degrees)
print(average_deg)
# compute r0 (TODO)
In [265]:
plt.figure()
pos = nx.spring_layout(graph)
nx.draw(graph, pos, node_color='b', node_size=10, with_labels=False)
plt.show()
Implementation of disease spreading¶
In [268]:
# set person class to each node
for i in range(len(graph)):
graph.node[i]["Data"] = Individual(i)
# create set of infecteds, initialize with 10 infecteds
init = np.random.choice(len(graph), 10, replace=False)
# time steps
t = 0
t_tot = 20
time_step = 365
dt = 1 / (time_step)
# constants of model
D = np.array([0.19,0.55]) * time_step
mu = 1e-3
PHI = [0.65, 0.50]
BETA = [0.59, 0.87]
# keep track of infecteds and steps
num_steps = int(t_tot / dt)
num_infected = np.zeros((2, num_steps))
num_res = np.zeros((2, num_steps))
p_inf = 0.8
# iterate over sextypes (MSM and HMW)
for sextype in range(2):
phi = PHI[sextype]
beta = BETA[sextype]
nu = (1 - phi)/D[sextype]
tau = phi/D[sextype]
beta = 1-(1-1/(average_deg*(1-p_inf)))**(nu)
num_res[sextype,:], num_infected[sextype,:], graph = network_model(beta, tau, nu, mu, init, num_steps, copy.deepcopy(graph), True)
Plotting¶
In [280]:
# plot number of infecteds
t = np.linspace(0,t_tot, num_steps)
plt.figure()
print("Mean number of infecteds:", np.mean(num_infected))
plt.plot(t, num_infected[0,:], label='all infecteds MSM')
plt.plot(t, num_res[0,:], label='resistant infecteds MSM')
plt.plot(t, num_infected[1,:], label='all infecteds HMW')
plt.xlabel("Time (years)")
plt.ylabel("number of people")
plt.plot(t, num_res[1,:], label='resistant infecteds HMW')
plt.legend()
plt.show()
In [281]:
plt.figure()
plt.plot(t,num_res[0,:]/num_infected[0,:], label='fraction of resistant strain MSM')
plt.plot(t,num_res[1,:]/num_infected[1,:], label='fraction of resistant strain HMW')
# plt.plot(num_res)
plt.xlabel("Time (years)")
plt.ylabel("Fraction of infected individuals")
plt.legend()
plt.show()
Animation¶
In [279]:
# set person class to each node
for i in range(len(graph)):
graph.node[i]["Data"] = Individual(i)
# create set of infecteds, initialize with 10 infecteds
init = np.random.choice(len(graph), 10, replace=False)
# time
t = 0
t_tot = 5
time_step = 365
dt = 1 / (time_step)
# constants of model
D = np.array([0.19,0.55]) * time_step
mu = 1e-3
PHI = [0.65, 0.50]
BETA = [0.59, 0.87]
# choose MSM or HMW to model
sextype = 0
phi = PHI[sextype]
beta = BETA[sextype]
nu = (1 - phi)/D[sextype]
tau = phi/D[sextype]
beta = 1-(1-1/(average_deg*(1-p_inf)))**(nu)
# initial state
_, _, graph = network_model(beta, tau, nu, mu, init, 0, graph, doInit=True, disable_progress=True)
# get values of disease status
values = []
for j in range(len(graph.nodes())):
if graph.node[j]["Data"].disease_status == 0:
values.append("Black")
elif graph.node[j]["Data"].disease_status == 1:
values.append("Blue")
elif graph.node[j]["Data"].disease_status == 2:
values.append("Red")
# create figure and set axes
fig = plt.figure(figsize = (10,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
im1, = ax2.plot([], [], c = 'b', label='all infecteds MSM')
im2, = ax2.plot([], [], c = 'r', label='resistant infecteds MSM')
ax2.legend()
ax2.set_xlim(0,t_tot)
ax2.set_ylim(0,num)
# initialize
pos = nx.spring_layout(graph)
plt.close()
ani_step = 365//5
frames = (t_tot*365)//ani_step
num_infected = np.zeros(ani_step*frames)
num_res = np.zeros(ani_step*frames)
# update function for animation
def update(i):
if i<1:
temp_res, temp_infected, _ = network_model(beta, tau, nu, mu, init, 1, graph, doInit=False, disable_progress=True)
else:
temp_res, temp_infected, _ = network_model(beta, tau, nu, mu, init, ani_step, graph, doInit=False, disable_progress=True)
num_infected[i*ani_step:(i+1)*ani_step] = temp_infected
num_res[i*ani_step:(i+1)*ani_step] = temp_res
t = np.linspace(0,i*ani_step/365, i*ani_step)
im1.set_data(t, num_infected[:i*ani_step])
im2.set_data(t, num_res[:i*ani_step])
ax1.clear()
# solution at step i
values = []
for j in range(len(graph.nodes())):
if graph.node[j]["Data"].disease_status == 0:
values.append("Black")
elif graph.node[j]["Data"].disease_status == 1:
values.append("Blue")
elif graph.node[j]["Data"].disease_status == 2:
values.append("Red")
nx.draw(graph, pos=pos, node_color=values, node_size=10, ax=ax1, clim = 2)
# Scale plot ax
ax1.set_title("Frame %d: "%(i+1), fontweight="bold")
ax1.set_xticks([])
ax1.set_yticks([])
# animate!
animation.FuncAnimation(fig, update, frames=frames, interval=100, repeat=True)
Out[279]:
Stochastic network model¶
The stochastic network model is based on an existing model by Kretzschmar, M. et al. (1996) [4]
The system consists of Persons defined by:
- Gender (male / female)
- Age (15 - 64)
- Sexual activity (high / low)
- Disease status (uninfected or (non-)resistant (a)symptomatic)
During each timestep, Partnerships are formed between two Persons. A Partnership is defined by:
- Type (casual or steady)
- The two Persons involved
The sexual contact network is therefore dynamic in time
Each timestep of the model consists of seven parts:
- Partnership formation
- Disease transmission
- Separation of partnerships
- Replacement
- Recovery
- Treatment
- Advance time
Initialize model¶
- Initialize 10,000 Persons of random age. 5% of Persons with age < 35 have high sexual activity
- Perform 2,000 iterations to obtain the steady state
- Infect 10% of the core group (high sexual activity) with non-resistant asymptomatic disease
The resistant disease enters the system through treatment; treatment will result in contracting the resistant disease with probability 0.01%
Time-evolution towards the steady state¶
In [231]:
fig
Time-evolution of the population¶
In [ ]:
fig
Time-evolution of infecteds over time¶
In [ ]:
fig
Fraction of resistant infecteds over all infecteds over time¶
In [ ]:
fig
Infectiousness of the disease ($R_0$)¶
Although the disease might seem quite infectious, simulations have shown $R_0=0.90 \pm 2.93$ (N=1,000, 95% CI). The spread of the disease is very dependent on where the disease is started. If a person with high sexual activity is infected, it will quickly spread throughout the entire core group whereas a person with low sexual activity may only affect a select few different partners
Animation¶
In [211]:
anim
References¶
[1] Fingerhuth, S. M., Bonhoeffer, S., Low, N., & Althaus, C. L. (2016). Antibiotic-resistant Neisseria gonorrhoeae spread faster with more treatment, not more sexual partners. PLoS pathogens, 12(5), e1005611.
[2] De, P., Singh, A. E., Wong, T., Yacoub, W., & Jolly, A. M. (2004). Sexual network analysis of a gonorrhoea outbreak. Sexually transmitted infections, 80(4), 280-285.
[3] Bansal, S., Grenfell, B. T., & Meyers, L. A. (2007). When individual behaviour matters: homogeneous and network models in epidemiology. Journal of the Royal Society Interface, 4(16), 879-891.
[4] Kretzschmar, M., van Duynhoven, Y. T., & Severijnen, A. J. (1996). Modeling prevention strategies for gonorrhea and chlamydia using stochastic network simulations. American Journal of Epidemiology, 144(3), 306-317.